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ABSTRACT 

We investigate the effects of magnetic fields and radiative protostellar feedback on the star 
formation process using self-gravitating radiation magnetohydrodynamical calculations. We 
present results from a series of calculations of the collapse of 50 M0 molecular clouds with 
various magnetic field strengths and with and without radiative transfer 

We find that both magnetic fields and radiation have a dramatic impact on star formation, 
though the two effects are in many ways complementary. Magnetic fields primarily provide 
support on large scales to low density gas, whereas radiation is found to strongly suppress 
small-scale fragmentation by increasing the temperature in the high-density material near the 
protostars. With strong magnetic fields and radiative feedback the net result is an inefficient 
star formation process with a star formation rate of < 10% per free-fall time that approaches 
the observed rate, although we have only been able to follow the calculations for 1/3 of a 
free-fall time beyond the onset of star formation. 
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1 INTRODUCTION 

Star formation is a remarkably inefficient process. This inefficiency 
in itself is a very good thing for the universe as a whole, since 
without it galaxies such as the Milky Way would very quickly ex- 
haust their supplies of gas by converting it into stars. Recent es- 
timates from the c2d Spitzer legacy survey of five nearby molec- 
ular clouds suggest that around 3 — 6% of the available gas in a 
molecul ar cloud is c onverted into stars in the local region of our 
Galaxy jEvans et alJi200 9). Previous observational results suggest 
simil a rly low efficiencies : e.g., ~ 1 — 6% in Taurus jEvans 
ll99lHOnishi etalJ|l998h : < 13% in the clouds in Chamaeleon 
jMizunoetal J 1999(1 though some dispers on in these results arises 
from the use of differing measures, whereas lEvans et afl use a uni- 
form definition of efficiency for all clouds. 

The source of such uniformly low efficiency is poorly under- 
stood, and it remains unclear as to what the "rate-limiting step" 
in star formation really is, since inefficiency is apparently present 
at a ll levels, from the formation of mole cular clouds in galax- 
ies toobbs et alJ l2008l : iLerov et alJ |200^ to the fact that only 
small, clustere d regions of molecular clouds with mass fractions of 
20 percent llLadal [1993 : Ijohnstone. Pi Francesco & Kirklliooi : 
iHatchell et aDl2005h participate in star formation, to the observa- 
tion that only a f raction of the mass in dense molecular clo ud cores 
ends up as stars jBenson & Mverslll989l : lAlves et al.l2007l) . Inspite 
of this, it is clear that a large part of the inefficiency lies within 
molecular clouds themselves. 



From a theoretical perspective, we have a very good idea of 
the basic ingredients of the star formation process - namely grav- 
ity, gas dynamics, turbulence, magnetic fields, radiative and me- 
chanical feedback, though their relative importance (particularly 
with respect t o magnetic fields and tu rbulence) remains vigourously 
debated (e.g. jCrutcher et alj I2OO9I : iMouschovias & TassisI I200A 
ICrutcher et al.l2008[) . By definition star formation involves the con- 
version of gas int o stars under self-gravity, the basics of which 
were elucidated bv lJeansI fT9o3) . The complication to the gas dy- 
namics is the highly turbulent (and supersonic) nature of molecular 
clouds and the wide range of length and time scales over which star 
formation takes place, presenting a formidable challenge for nu- 
merical simulations even before considering other relevant physics. 
Neverth eless, simulations including just self-gravity and hydrody- 
namic s I Klessen. Burkert & Bate! 1 19981: iBate. Boimell & Broirmi 



Bate 



2002allbll2003l : lBonnell. Bate & Vinell2003l : lBate & Bonnellll2005l: 



2005h have been surprisingly successful in predicting many 



properties of clustered star formation, including the initial mass 
function (though with an overproduction of brown dwarfs), mul- 
tiplicity as a function of primary mass, the frequency of very low 
mass binaries, general trends for the separation and mass ratio dis- 
tributions of binaries and the relative orbital orientations of triple 
systems (Bate 2009a). 

However, the efficiency of star formation in these calculations 
would be J> 50% were the simulations left to run, since in the ab- 
sence of stellar feedback there can be nothing to prevent all of the 
(bound) gas from eventually accreting onto the stars. While tidal 
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forces from the underlying galactic potential may be a contribut- 
ing factor in some cases (Ballesteros-Paredes et al. 2009), many 
suggestions in the literature are related to the turbulence present 
in the cloud. Star formation in initially unbound clouds, where 
star formation wil l nevertheless occur in the presence of a turbu- 
lent velocity field dClark & Bonnell" 2004;' Clark et al."200 ^, is one 
extreme. But sources of turbulent driving (e.g. Matzner & McKed 
hood iKrumholz. Matzner & McKeell2006l : iNakamura & Lill2007h . 
whilst not strictly changing the overall efficiency in a bound cloud, 
can dramatically alter the fraction of a cloud which is unstable to 
gravitational collapse in a given dynamical time, decreasing the ef- 
ficiency per free-fall time (e.g. |Padoanlll995l : iKlessen et al]|200d : 
iKrumholz & McKeel2005h . This occurs naturally in turbulent mod- 
els because there is a spectrum of density fluctuations, of which 
only a small fraction is in sub-regions dense enough to be Jeans- 
unstable. Star formation is thus made inefficient, in a per free-fall 
time sense, because turbulence produces in a given dynamical time 
a range of clu mps (or "cores"), on ly some of which are bound and 
will collapse i Klessen et al. There is still a difficulty, how- 

ever, which is that if this is to work for dense, globally bound re- 
gions, one has to keep driving the turbulent motions (either from 
outside or within) otherwise the turbulence will quickly decay and 
the gas from the whole dense globally bound clump will eventually 
be used up in star formation. 

Magnetic fiel ds have long been recognised as a key ingredient 
in sta r formation jMestel & Spitzeilll956l : IShu et alJ[T987l : iMestell 
1 19991) . given that observations robustly measure fields at sufficient 
strengths that they are close to preventing star formation altogether 
dMac Low & KlessenI liooj) . and robustly in th e regime where 
magnetic pressure is dominant over gas pr essure ( ICrutchei|[T99^ : 
iBourke et al.l200ll : lHeiles & Trolandll2005l) . The importance of the 
latter point is easily overlooked and implies that, even if magnetic 
fields do not prevent global gravitational collapse in a molecular 
cloud, they can neve rtheless act as the dominant source of pressure 
( |Price&Batell2008ll . supporting la rge fractions of the clo ud and 
perhaps regulating star formation jNakamura & Lil l2005h . Mag- 
netic fields, may also have important effects on the statistics of 
molecular cloud turbulence, even in the regime whe re the Alfven 
speed is small compared to the turbulent velocities jPadoan et al.l 
|2007|) . Importantly, magnetic fields are not usually included in de- 
terminations of whether or not a molecular cloud core is "gravita- 
tionally bound". 

Radiation presents a complementary method for regulating 
star formation and the need for star formation simulations to 
incorporate the ef fects of radiati ve transfer has a lso long been 
understood (e . g. iLarsor] Il969l: Is iack & Boden heimej 1 19751 : 
IBoss & Mvliilil 1 19921 : iMasunaga & InutsukiJ \2o6dtl . Radiative 
feedback affects star formation as soon as the gas becomes 
optically thick, setting the "opacity Unfit" beyond which frag- 
mentation can proceed no further jLow & Lvnden-Belll 1 19761 : 
^ees 1976). From thereon the newborn protostar can continue 
to radiate into the surrounding gas, increasing t he Jeans mass 
and thus inhibiting further star f ormation ( Whitehouse & Batd 
I2OO6I : lKrumhoi3l2006l: lBatell2009bl) . In the case of massive stars 
radi ation may be sufficient to halt accretion from the cloud ( iKahnI 
1 1974 : Wolfire & Cassinelli 1987) although variou s non-spherical 
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rithms f or magnetohydrodynamics (MHD) jPrice & Monaghari 
l2004allbL [200^: iPrice &Batdl2007h and radiative transfer in the 
flux-limited diffusion approximation dWhi tehouse & Batd l2004l : 
IWhitehouse. Bate & Monaghad l2005l : [whitehouse & Batd I2OO6I) 
in the context of the smoothed particle hydrodynamics (SPH) 
method has nevertheless paved the way for such effects to be in- 
cluded. Thus we have recently been abl e to incorporate, th ough 
separately, the effects of magnetic fields (I Price & Bat^bOO Sl) and 
radiative feedback ( Bate 2009b) into simulations of star cluster for- 
mation. In this paper we study, for the first time, the combined ef- 
fects of both. 

We thus present a 'recipe' for inefficient star formation in even 
relatively dense molecular clouds. The ingredients (i.e., the equa- 
tions of self-gravitating radiation MHD and our numerical formu- 
lation of them) are presented in Section|2] The initial conditions for 
our simulations are discussed in Section [3] We present our results 
in SectionUand discuss their wider implications in Section|5] 



2 NUMERICAL METHOD 

In this paper we solve the equations of self-gravitating radiation 
magnetohydrodynamics, using a two-temperature flux-limited dif- 
fusion scheme for the radiation, coupled with the equations of 
ideal MHD (that is, assuming infinite conductivity and without 
considering ambipolar (ion-neutral) diffusion or the Hall effect). 
Whilst we have previously published st ar cluster formation cumu- 
lations using separately either the MHD ice & Bate"2007','200d) 
or usin g the flux-limited diffusion (Whitehouse & Bate 2006; Batg 
l2009bb schemes, this is the first time which we have combined the 
t wo. Thus, whilst the MH D formulation is ide ntical to that used 
in IPrice & Batd ( l2007l) and lPrice & Batd j2008h and the radiation 
scheme is based on that used in l Batd ( l2009b ir some minor changes 
have been made to the radiation terms in order to combine them 
with the MHD part of the code. 



2.1 Equations of Radiation Magnetohydrodynamics 

The equations of self-gravitating radiation MHD are solved in the 
form 
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However, neither magnetic fields nor radiation are easy 
to incorporate into three-dimensional numerical simulations 



where p is the density, v is the velocity, P is the hydrodynamic 
pressure, B is the magnetic field, u is the specific thermal energy 
of the gas, $ is the gravitational potential, ^ and Prad are the 
frequency-integrated specific radiation energy and radiation pres- 
sure tensor respectively; a, c, x, i<- and c„ are the radiation constant, 
the speed of light, the total and absorption opacities and the ratio of 
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specific heats respectively and F is thie radiative flux, which in the 
flux-limited diffusion approximation is given by: 

F = — V(pC), (8) 
up 

where A is the dimensionless flux limiter which is designed to en- 
sure thatjhe radiation propagates no faster than the speed of light 
(see lWhitehouse et alj|2005l for details). The above expression for 
F means that the first term in equation l[4j becomes a diffusion term 
for the radiation energy (hence "flux-limited diffusion"). 

Equation (TJ is an exact solution to the continuity equation 
which is re presented in SPH form by the density summation (see 
|Pricell2008l for the difference between integral and differential for- 
mulations in an SPH context). Equation (O is the equation of mo- 
tion for the gas which contains force terms from the hydrodynamic 
(VP), magnetohydrodynamic V (5^^ ^ ■'^■'^)] '"^'^ radia- 
tion (F, ie. V(p5)) pressure gradients and from the gradient in the 
gravitational potential (V^"). Equations l|3j and l|4j are the energy 
equations for the gas and radiation respectively. Equation l|5j is an 
expression of the magnetic field in terms of the Euler (or Clebsch) 
potentials as and Pe which maintains the divergence constraint 
(V ■ B = 0) by construction and for which the induction equa- 
tion for the mag netic field takes the particularly sim ple form given 
by equation l|6j JSternll 19701 : iRosswog & Pricell2007t) . It should be 
noted that use of the Euler potentials approach also introduces lim- 
itations on the topology of fields that can evolve d uring the cal- 
culatio n. Whilst thes e are discussed in more detail in lPrice & Bate! 
l l2008h and lPrice & B ate ( 2007), the main physical process not cap- 
tured is the winding up of magnetic fields on smaller scales, since 
the Euler potentials rely on a well defined mapping from the initial 
particle positions to those at a later time. This means that, whilst we 
are able to study the influence of magnetic fields on the large scale 
structure of the cloud, field growth on smaller scales is not well cap- 
tured. On the other hand the ideal MHD approximation also breaks 
down at these scales, so an improved formulation would also need 
to correctly account for non-ideal MHD effects such as resistivity 
and ambipolar diffusion. Finally, Poisson's equation (Equation |7) 
is solved in order to determine the gravitational force. 

The equation set is closed by equations of state for the gas and 
the radiation field. For the gas, the equation of state is given by the 
ideal gas law 

P=^, (9) 

where T is the gas temperature, TZ is the gas constant and jj, is the 
mean molecular weight. The equation of state takes into account 
the translational, rotational and vibrational degrees of freedom of 
molecular hydrogen (assuming a 3: 1 mix of ortho- and para- hydro- 
gen th at remains fixed throughout the calculations; see lBolev et al.l 
l2007f ). It also includes the dissociation of molecular hydrogen and 
the ionisations of hydrogen and helium (which are assumed to have 
mass fractions of X — 0.7 and Y — 0.28 respectively). The con- 
tributions of metals to the equation of state is neglected. 

For the radiation, the equation of state is given by the Edding- 
ton approximation 

Prad = fpC, (10) 

where f is the Eddington tensor which has both an isotropic term 
and an anis otropic term related t o the gradient in radiation energy 
density (see Whitehouse & Bate 2006 for details). 

For comparison with previous results, we have also performed 
a set of calculations without radiative transfer, but which use a 



barotropic equation of state for the gas (i.e. replacing equations [3] 
|4]and|9} of the form 

P = Kp''. (II) 

where the polytropic exponent 7 is given by 

7 = 1, p< IQ-^-'g cm-^ 
7 = 7/5, p > cm^^ (12) 

The simulations using the baro tropic equation of s tate are thus 
identical to those performed by IPrice & Bate! ^20081) except for a 
factor-of-ten decrease in the sink particle radii (see below) and also 
the MPI-parallelisation of the tree-code, both of which change the 
overall fragmentation pattern slightly due to the chaotic nature of 
star formation. 

2.2 Numerical method 

We solve equations l|T](-([7} using the Smoothed Particle Hydro- 
dynamics (SPH) meth od (for reviews see Monaghan 1992; Pric^ 
|2004 lMonagharj|2005 ) . The SPH formulations of various parts of 
these equations, as currently implemented in our code, have been 
separately described and tested in a number of papers (mostly in- 
volving the authors), as summarised below. 

2.2.1 SPH formulation 

The self-gravitating part of the algorithm (i.e. equation Q and the 
gravitational force term in equation {2)) is ident ical to the energy- 
conser ving formulation described and tested in iPrice & MonaghanI 
( l2007h . The gravitational force is softened using a softening length 
that is equal to the SPH smoothing length and formulated such that 
taking the Laplacian of the gravitational potential results in pre- 
cisely the right hand side of Poisson's equation (O with the density 
p equal to that calculated in the hydrodynamic s via the SPH sum- 
matio n (ie. the SPH expression of equation[l]l jPrice & Monaghai] 
'2007"). Furthermore - despite the softening length being a variable 
function of position - momentum, energy and angular momentum 
are conserved exactly using this formalism. However, in practice, 
a nearest-neighbour binary tree algorithm is used to efficiently cal- 
culate the long-range part of the gravitational force (and also return 
the list of SPH neighbours), which does not conserve momentum, 
angular momentum or energy exactly. The tree code formed the 
original core of the SPH code and remains essentially as originally 
implemented bv lBenz et alj ( Il990h . 

All of the evolution equations are integrated using a second 
order Runge-Kutta-Fehlberg method with a timestep related to the 
convergence of each variable with timestep (that is, determined by 
comparing the error using half of the current timestep to that over 
the full timestep). Given the rich array of physics in our current 
calculations, we have found this to be substantially more accurate 
than a standard leapfrog method where the timestep is based only 
on stability considerations (rather tha n accuracy). Individual par- 
ticle timesteps were implemented by iBatd ( Il995h in order to ef- 
ficiently follow calculations where the timestep is constrained by 
only a small fraction of particles in a simulation. 

The hydrodynamics and MHD parts of the code (i.e., 
the numerical formulation of equations l[T}-l|3j apart from 
the radiation terms and equations ((S}-©) are based on the 
smoot hed particle magnetohydrodynamics algorithm developed by 
IPrice & Monaghan ( 2004a, b, 2005) and ap plied to star forma- 
tion using the Euler potentials formulation bv |Price&Bat3 ( l200l 
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I2OO8I ) (see alsplRos swog & Pricj20o"7h . Special attention has been 
paid to the formulation of terms relating to the gradient of smooth- 
ing length jPrice & MonaghanI l2004bL l2007h which ensures that 
total energy and entropy are conserved exactly by the hydrody- 
namic parts of the equations. Energy and momentum conserva- 
tion is not maintained exactly for the magnetic parts of the equa- 
tions in order to avoid the well-known instability relating to ex- 
actly mo mentum-conserving form ulations of the SPMHD force 
term fsee lPrice & MonaghanI |200^ . Furthermore, using the Euler 
potentials, the force equation is not directly derived from the nu- 
merical form of the induction equations leading to a very small 
error in energy conservation. In practice, however, these errors are 
much smaller than those introduced by the treecode for the gravita- 
tional force and individual particle timesteps. 

Dissipative terms corresponding to artificial viscosity and ar- 
tificial resitivity are added in order to capture shocks and mag- 
ne tic reconnection, re s pective ly. These are applied as described 
Price & Monai and for th e Eule r potentials by 



IPrice & Bate! ( l2007h and lRosswog & Pried ( l2007h . For the calcu- 
lations with radiative transfer, the energy associated with this dissi- 
pation is added to the thermal energy, though obviously such energy 
is discarded for the calculations employing a barotropic equation of 
state and thus also during the initial period {t < It//) for all the 
calculations during which the barotropic equation of state has been 
used (see below). We are, therefore, not able to realistically assess 
the effect of any heating that may arise due to magnetic reconnec- 
tion since any heat created by reconnection in the early phases is 
lost and at later times, whilst the energy is captured, the field struc- 
ture is effectively lost because of the limitations to the Euler poten- 
tials approach on smaller scales (see above). Thus, the contribution 
of magnetic dissipation to heating in the present calculations is very 
small. 

The radiative transfer parts of equations l|2j and ([3} and 
the radiative energy equation ^ are solved implici tly using the 
for mulation develope d bv jWhitehouse & Bate! ( |2004) . accelerated 
by Whitehouse et al. 1 200^ a nd ap p lied to star formation by 
IWhit ehouse & Bate (2006^ and Bate ('2009b'). We use the same 
opacities as Whitehouse & Bate (2006). In order to combine the 
radiat ive transfer parts of the code (developed bv lWhitehouse et al.l 
I2OO5I and based on a traditional "number of neighbours" approach 
to variable smoothing lengths in SPH) with the MHD (developed 
bv lPrice & Monaghanll2004bl . which is formulated taking account 
of variable smoothing length gradient terms), minor modifications 
have been made to the manner in which the radiative transfer equa- 
tions are expressed in SPH form. The main change is that the ra- 
diation diffusion term in equation lO is calculated using an aver- 
age of the kern els rather than an average smoothing length as in 
I Whitehouse & Bate (2006|), so that the diffusion term in the energy 
equation becomes 



where the average of the SPH kernel gradients is 
1 



vw,j = -[vw,Ah) + vw,,ih,)], 



and 



(14) 



(15) 



is in turn a function of smoothing length via the SPH density sum- 
mation. A solution to the densit y summation must therefor e be ob- 
tained iteratively as described in lPrice & MonaghanI ( l2007h . 



2.3 Sink particles 

Sink p articles were introduced into SPH bv lBate. Bonnell & Pried 



Sink p a 
( 11993) 



in order to follow star formation calculations beyond the 
formation of the first star. 

For the calculations presented here that do not include radia- 
tive transfer, the cri terion used fo r sink p article creation is identical 
to that described in lPrice & Batd ( l2008 h and we therefore refer the 
reader to that paper for deta ils. The major difference between the 
calculations of lPrice & Batd and the similar calculations in this pa- 
per is that we have used sink parti cles with an accretio n radius of 
only 0.5 AU, compared to 5 AU in IPrice & Batd i l2008h (and sim- 
ilarly in BBB03). This adds considerable computational expense 
to the calculations because the closest gas orbits around the sinks 
(and the smallest length scale ~ h) are reduced by a factor of 10, 
resulting in an increase in the maximum acceleration by a factor 
of — 100 and therefore a decrease in the minimum timestep 
in the calculations At cx ^/h/\a\ by a factor of ~30. Whilst such 
expense is unnecessary when using a barotropic equation of state 
(in Section [4| we comp are our results to previous results obtained 
bv |Price&Batdl2008l using 5 AU sink radii and find essentially 
no difference to the fragmentation), it is important for the radiative 
transfer calculations. 

With the exception of the addition of magnetic fields, the ra- 
diative transfer calcul ations presen ted here are very similar to those 
recently published bv lBatd ( l2009bl) . In both, the gas is followed be- 
yond the first hydrostatic core phase ( lLarsoiilll969h and the onset 
of molecular hydrogen dissociation (T « 2000 K). Sink particles 
are inserted during the second collapse pha se, ju s t before a stellar 
core would be formed in the calculations. iBatd ( l2009bh inserted 
sink particles at a density of 10~^ g cm"'' while in the calculations 
presented here they are inserted slightly earlier at 10~^ g cm~^. 
In terms of the real star formation process, t his is just a co uple of 
weeks before the stellar core is formed. As in lBatd ( l2009bl) . no ra- 
diative feedback is provided by the sink particle. The radiative feed- 
back provided by the protostars is limited to the radiation emitted 
from the gas as it falls into the sink particles. Thus, it is important to 
make the sink particle accretion radii as sm all as is com putationally 
practical (0.5 AU in both the calculations o flBatel2009bl and the cal- 
culations presented here). As noted bv iBatd ( l2009bh . because not 
all of the protostellar luminosity is fed back into the calculations, 
the effects of radiative feedback seen in the calculations presented 
here must be viewed as a lower limit. 

An estimate of the energy input that we are missing from ac- 
cretion within the sink particle radius can be made by comparing 
the accretion luminosity from within this region to that available 
from infinity. The accretion energy expected from within our sink 
particle radius of 0.5 AU is given by 



(16) 



Note that in the variable smoothing length formulation of SPH, the 
smoothing length h is an analytic function of the density p, which 



This may be compared to the accretion luminosity that is captured 
in our calculations by accretion to the sink radius from infinity, 
which is given by Lace ~ GAIM /{0.5 AU). Therefore, if we as- 
sume a protostellar radius of ~ SRsun there is potentially a factor 
of up to ~ 30 in further energy input that is missing from the cur- 
rent calculations. Bate (2009b) investigated the effect of this miss- 
ing radiation on his similar calculations that did not include mag- 
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netic fields by repeating a calculation with a larger accretion radius 
of 5 AU (i.e. reducing the accretion luminosity by a further order 
of magnitude). He found only a small difference in the amount of 
fragmentation that occurred between the 0.5 AU and 5 AU calcula- 
tions because even the heating present in the calculation with 5-AU 
accretion radii was enough to inhibit fragmentation near to existing 
protostars. Thus, while we again emphasize that the radiative feed- 
back incorporated into the current calculations is only a lower limit, 
we believe that using accretion radii of 0.5 AU captures the essence 
of the effects of radiative feedback, at least in terms of fragmenta- 
tion. This is in stark contrast to the situation encountered using a 
barotropic equation of state. 



3 INITIAL CONDITIONS 

The initia l conditions for the s imulations are identical to those pre- 
sented by^ ^^^^^^3 HqO^) and similar to the original calcula- 
tio n of | Bate et al ]( l2003h and the first of the calculations performed 
bv lBatdbOOm We briefly recap the initial conditions below. 



3.1 Density, temperature and velocity field 

We set up an initially uniform, spherical cloud with a diameter of 
0.375 pc (77,400 AU) that contains a total of 50 Mq of molecular 
gas, giving an initial density of po = 1.2 x 10~^^g cm~"^ (nn^ = 



3 X 10'*) and a global free-fall time of ts = y^3n/ {32poG) = 
1.90 X 10^ yrs. 

The cloud is constructed using 3.5 million SPH particles (de- 
ter mined by the resolutio n re quirement for reso lving the Jeans mass 
bv lBate & BurkertlliggTl . see Bate et aDl2003h placed in a uniform 
random distribution cropped to the cloud radius (i.e. no particles 
are placed exterior to the cloud). This results in a significant expan- 
sion of the outer layers as the calculation proceeds (equivalent to 
the assumption of open boundary conditions in a grid-based sim- 
ulation). The initial sound speed was set to 1.84 x 10* cm s^^, 
corresponding to a temperature of 10 K given the mean molecular 
weight of p ^ 4.0/(2 x 0.7 + 0.28) = 2.38 amu. The resultant 
ratio of thermal to gravitational energy was Ograv = 0.074. 

A supersonic 'turbulent' velocity field with power s pectrum 
P(k) oc A;"* (i.e. consistent with Larson's scaling relations. lLarsoij 
198 1J) was impos ed upon the initially uniform density cloud as in 



Bate et al ] ( l2003h . with the initial velocity field normalised such 
that the kinetic energy is initially equal to the gravitational potential 
energy of the cloud. This gives an initial root mean square (RMS) 
Mach number of 6.4 and an initial RMS velocity of 1.17 x 10''' 
cm/s. 

The computational challenge of star formation is well illus- 
trated by the fact that during the calculations, we find that the dens- 
est regions can contain particles moving on a timestep up to 2^' 
times smaller than the largest timestep bin (which is constrained 
by the time between output dumps), so that the shortest timestep 
is around 1.5 hours compared to a total evolution time of several 
hundred thousand years. 



3.2 Magnetic fields 

We quantify the relative strength of the magnetic field in terms of 
the mass-to-flux ratio (i\//<l>) of the cloud, compared to the critical 
value for the onset of collapse in a spherical cloud given by (e.g. 



iMou schovias & Spitzer|[l97i : lMestellll999l : | Mac Low & Klessei] 
,2004) 



2ci 
3 



TvGpo 



(17) 



where G and fio are the gravitational constant and the permeability 
of free s pace respectively and ci is a c onstant determined numeri- 
cally bv lMouschovias & Spitzerl l ll976h to be ci ^ 0.53. 

In this paper, we have performed calculations starting with an 
initially uniform magnetic field with mass-to-flux ratios in units 
of the critical value of i\f/$ = oo (i.e. no magnetic field), 10, 5 
and 3. All of our calculations are 'supercritical' (that is, unstable to 
collapse) as under our assumption of ideal MHD (i.e. no ambipolar 
diffusion or resistivity), subcritical clouds would not (and do not) 
collapse. 

The corresponding physical field strength for a given mass-to- 
flux ratio and cloud dimensions is 



Bo = 194 pG ( ^ 



M 



50 Me 



R 



0.188 pc 



(18) 



where Af/$ is the mass to flux ratio in units of the critical value. 
Thus, a simulation with a critical mass-to-flux ratio would have 
-Bo ~ 194/^G and for the calculations with mass-to-flux ratios of 
oo, 10, 5 and 3 the corresponding field strengths are given by Bo = 
0, 19, 39 and 65/iG, respectively. 

The magnetic field may also be parametrised in terms of the 
plasma the ratio of gas to magnetic pressure, according to 



P = 0.028 ( f 



Vl84ms-i/ 



M 



50 Mr. 



R 



0.188 pc 



(19) 

The simulations presented here thus have initial /3's of oo, 2.8, 0.7 
and 0.25 respectively. Note that the magnetic pressure is dominant 
over gas pressure in the cloud for mass-to-flux ratios < 6 which 
is the case for the two strongest-field calculations. Indeed, as in 
IPrice & Bat3 ( I2OO8I) . we find that these two calculations show far 
more significant differences compared to the weaker field and hy- 
drodynamic calculations. 

Finally, the Alfven speed in the initial cloud can be computed 
using 



VA = 1.6 X 10^ cms ^ 



M 



50 Mr; 



R 



0.188 pc 



(20) 

giving VA = 0, 1.6 x lO'', 3.1 x 10'' and 5.2 x 10'' cm s~^ for the 
calculations in this paper. Thus, the initial turbulent motions in the 
cloud are super- Alfvenic in all cases with Alfvenic Mach numbers 
of 00, 7.3, 3.8 and 2.3, respectively. 



4 RESULTS 

We have computed a total of eight calculations, that is, for four 
different mass-to-flux ratios, both with and without radiative trans- 
fer (where "without" means that we use the barotropic equation of 
state given by equation l ll2t instead). 

The evolution of the simulations can be divided into two 
stages: i) the initial collapse of the cloud (i.e. up to ~ 1 free-fall 
time) during which the cloud is optically thin, essentially isother- 
mal, and the dynamics and large scale structure are primarily con- 
trolled by the interaction of turbulence and magnetic fields; ii) the 
subsequent evolution of the cloud after the formation of the first 
star (i.e. > 1 free-fall time), where the cloud has optically thick 
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Figure 2. A comparison of the fragmentation that has occurred in the eight different calculations at 1.2 free fall times (iff). The two left-hand columns show 
column density for each calculation, with magnetic field strength increasing from top to bottom (as indicated by the mass-to-flux ratio in units of the critical 
value for collapse, where M/<i> = oo corresponds to hydrodynamics) using either a barotropic equation of state (first column) or with radiative transfer 
(second column), as indicated. A strong decrease in protostar formation with increasing magnetic field strength may be observed (comparing rows from top to 
bottom). The radiative feedback from the protostars is illustrated by plots of the mass-weighted temperature (J pT dz/ f p dz), shown in the corresponding 
right-hand panels. The effect of the radiation heating the gas in the vicinity of the protostars (fourth column) can be seen to be poorly captured by the barotropic 
equation of state approximation (third column) and leads to a dramatic suppression of small-scale fragmentation (comparing the first and second columns). 

regions embedded in the wider (optically thin) large scale struc- 4.1 Large-scale cloud structure 
ture, and where the small-scale fragmentation is regulated by the 

radiative feedback from existing protostars on the gas. During the first phase, the radiative transfer has little effect on the 

overall dynamics compared to the use of a barotropic equation of 
state because the cloud is optically thin, radiation can escape eas- 
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ily and there are no significant sources of radiation. Typical tem- 
perature variations are of the order AT/T ^ 10%. However, it 
is computationally very expensive to compute the evolution of the 
cloud with radiative transfer in the optically thin regime. We have 
therefore computed only one set of calculations of the global cloud 
structure during the first (isothermal) phase, the results of which are 
shown at one free-fall time in Figure[T](with magnetic field strength 
increasing from left to right, as indicated). These are essentially the 
same as those presented bv lPrice&B"at3 ( l2008h are we therefore 

discuss them only briefly here. 

Figure[T] as in |Price&Batel ll2008^ reveals the dramatic influ- 
ence the global magnetic field has on the large scale structure of 
the cloud, even though the field is much too weak to prevent global 
gravitational collapse. In particular, for the two strong magnetic 
field calculations (mass-to-flux ratios of 5 and 3 shown in the two 
right hand panels) large-scale voids are visible in the cloud where 
material has slipped down the field lines to leave behind evacu- 
ated but magnetically-pressurised voids. The se magnetic-p r essure 
supported voids were discussed in detail in ' Price & Batj j2008l) 
(see also Price et al. 2008) and appear in the regime where /3 < 1 
(i.e. where the magnetic pressure is dominant over the gas pres- 
sure). This regime is particularly interesting given that almost all 
magnetic fie ld strength mea s urements in molecular clouds indicate 
that < 1 (Crutcher 1999: 'B ourke et al.ll200ll : lHeiles & Troland 
|2004 iHeile s & Crutcher 2005). Also visible during the initial ex- 
pansion phase is a 'stripy' structure in the column density maps 
which is aligned with the large scale magnetic field lines. This is 
a consequence of the anisot ropy of turbulent motions in the pres- 
ence of a magnetic field (e.g. lGoldreich & Sridha r 1995) and, while 
not so obvious in Fi gurelH was discussed and clearly illustrated by 
IPrice & Batdtosl) . 

4.2 Fragmentation 

From one free-fall time (tg), the cloud structures shown in Figure[T] 
were evolved both with and without radiative transfer (i.e. using 
the barotropic equation of state in the former case and the full flux- 
limited diffusion equations in the latter). The simulations were run 
from this point to between 1.25 and 1.54 ta (2.93 x 10^ yrs) de- 
pending on the computational expense (the calculations slow down 
significantly once star formation initiates and the more protostars 
are formed, the slower the calculations become). The barotropic 
calculations with mass-to-flux ratios of A//$ = oo, 10, 5 and 3 
begin forming stars ai t ~ 1.07, 1.03, 1.10 and 1.19tfi, respec- 
tively, with the star formation in the radiative transfer conterparts 
typically being delayed by ~ 0.01 tg. A close-up of the fragmenta- 
tion in all eight simulations is shown in Figure |2l showing column 
density (left-hand panels) and mass-weighted temperature (right- 
hand panels) at 1.20 f ft , after star formation has begun in all eight 
clouds. The sink particles are shown as white filled circles. 

The left-hand (column density) panels of Figure |2] dramati- 
cally illustrates two main effects. The first is an overall decrease 
in star formation rate with increasing magnetic field strength (rows 
from top to bottom are in order of weakest to strongest magnetic 
field). This is a result of the influence of the global magnetic field 
on the large scale cloud structure, as already evident in Figure [T] 
In particular, for the stronger field calculations (bottom two rows 
of Figure IJl and the rightmost two panels of Figure [T}, large parts 
of the cloud are supported against collapse by the magnetic field 
resulting in fewer collapsing sub-regions (or 'cores'). For example, 
where the hydrodynamic calculation (top row of Figure |2]l has col- 
lections of protostars separated by a couple of thousand AU, only 



one collapsing region is evident in the strongest field (M/<1? = 3) 
case (bottom row), which shows no sub-fragmentation either with 
or without radiative transfer. The effect of the magnetic field in 
slowing the infall from the global cloud is further quantified in Fig- 
ure[3]and discussed in Section l43] below. 

The second effect visible in Figure |2] is the dramatic suppres- 
sion of small-scale fragmentation by the radiative feedback. This 
is especially obvious in the hydrodynamic/weak field calculations 
where the calculations using a barotropic equation of state have 
fragmented into multiple low-mass objects which interact violently, 
causing ejections of very low mass objects from multiple systems. 
This small-scale fragmentation occurs primarily in the massive pro- 
tostellar discs. By contrast, in the calculations which include radia- 
tive transfer, all of the subsequent disc fragmentation is suppressed 
by the radiation from the existing protostar(s). A good example is 
found in the lower-right of the M/$ = 10 panels of FigurelJ] with 
radiative feedback a single object with a disc is formed, while with- 
out radiative feedback this disc fragments into three objects, one 
of which is ejected. The radiative feedback, in effect sets a min- 
imum distance between protostars by substantially increasing the 
temperature and therefo re the Jeans length in the gas immediately 
surrounding a protostar ( lBatel2009bh . 

The differences between computing the radiative transfer and 
using the barotropic equation of state approximation are best illus- 
trated by plotting the temperature, given in the right-hand panels 
of Figure|2] Each panel shows the integrated temperature map (i.e. 
J pT dz/ J p dz) for the corresponding column-density panel in 
Figure|2] We have not plotted the sink particles on these panels so 
that the temperature distribution very close to the protostars can 
be seen for the barotropic calculations. For the barotropic equa- 
tion of state (centre-right panels), the temperature is simply related 
to the density, leading to very point-like sources of energy con- 
centrated around the protostars themselves. In the radiative trans- 
fer calculations (right-most panels), the radiation emitted from the 
high-density optically-thick gas near the protostars heats a much 
larger surrounding region to temperatures > 30 K, effectively shut- 
ting off any further fragmentation in this material (as evident in the 
left-hand column-density panels of Figure[2ll. 



4.3 Star formation rate 

The effects of both magnetic fields and the radiative feedback on 
the star formation rate are quantified in Figure [3] which shows the 
total mass in protostars (that is, the total mass of all sink parti- 
cles in a simulation) as a function of time. After t = 1.2 tg, 
the eight simulations form a strict sequence of progressively de- 
creasing star formation rate in the order; hydrodynamic, barotropic; 
hydrodynamic, RT; M/$ = 10, barotropic, M/$ = 10, RT; 
M /$ = 5, barotropic, M/$ = 5, RT; M/$ = 3, barotropic, 
Af/<1> — 3, RT, i.e. with magnetic fields as the primary effect and 
radiative feedback secondary. The rate at which gas is converted 
into stars decreases due to the influence of both magnetic fields 
and radiative feedback, though more strongly with the former. For 
example, at t — 1.2 tg, the hydrodynamic calculation contains 

2.2 Mq in stars using a barotropic equation of state, compared to 

1.3 Mq with radiative feedback, both of which are higher than the 
[0.83 Mq , 0.44 Mq] formed by the weak field M /$ = 10 calcula- 
tion at the same time [without,with] radiative transfer. These num- 
bers decrease further to [0.42 Mq,0.33 Mq] for the M/$ = 5 
simulation and further still to [0.055 Mq, 0.0056 Mq] for the 
strongest magnetic field case (M/"I> = 3). This general trend is 
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Time [Initial free-fall times] 

Figure 3. Total mass in stars (sink particles) as a function of time, show- 
ing all eight calculations with (thick lines) and without (thin lines) radia- 
tive transfer at four different magnetic field strengths: hydrodynamic (solid 
black lines), Af/<3? = 10 (dotted red lines), M/$ = 5 (dashed blue lines), 
and M/^ = 3 (dot-dashed magenta lines). The star formation rate de- 
creases with increasing magnetic field strength and with the addition of ra- 
diative feedback. Note how the two curves for each magnetic field strength 
track each other for some time before diverging, indicating that radiative 
feedback only plays a role in suppressing subsequent fragmentation rather 
than changing the initial pattern of star formation. 



continued as far as we have been able to run the calculations in 
each case (Figure[3](. 

The fact that the radiative feedback influences subsequent star 
formation rather than the initial fragmentation is also evident from 
Figure [3] In particular, the two curves corresponding to the same 
magnetic field strength but with and without radiative transfer in 
each case track each other closely after first sink formation, before 
diverging at later times. Taking the Af/<1> — 3 case as an example 
(i.e. the lower two curves in Figure[3}, and comparing the time evo- 
lution in Figure[3]to the fragmentation sequence shown in Figure[5] 
it may be observed that the two curves diverge when secondary disc 
fragmentation occurs in the barotropic calculation (t ~ 1.27 tft), 
leading to a burst of star formation (and subsequent ejection of low- 
mass objects from the multiple system). In the radiative transfer 
case, the disc does not fragment but instead continues to slowly 
accrete onto the existing protostar. 

Magnetic fields and radiation are also found to affect dif- 
ferent densities in the cloud. Figure |4] shows the mass above a 
given density threshold in the cloud as a function of time for 
three different density thresholds, p > 10~^^g cm~"^ (bottom 
panel), p > lO^^^g cm^'^ (middle panel), and p > 10~^^g cm~"^ 
(top panel), where solid lines correspond to calculations using a 
barotropic EOS and dashed lines refer to calculations using radia- 
tive transfer and, as in Figure [3] the lines form a sequence from 
top to bottom with increasing magnetic field strength. At a den- 
sity threshold of lO'^^g cm""^ (bottom panel), whilst there is a 
strong decrease in the mass collapsing to higher densities with in- 
creasing magnetic field strength, there is almost no difference be- 
tween the barotropic simulations and those with full radiative trans- 
fer (i.e. comparing the solid and dashed lines), indicating that radia- 
tive feedback plays very little role at these densities. At a threshold 
of p > lO'^^g cm~"^ (middle panel) the results are similar (al- 




0.8 1 1.2 1.4 1.6 



Time [Initial free-fall times] 

Figure 4. The total mass above certain density thresholds in each col- 
lapsing cloud as a function of time. From top to bottom the panels show 
M{p > 10~^^g cm"'') (approximately two orders of magnitude denser 
than the original cloud density), M{p > lO^'^^g cm"''), and M(p > 
lO^^^g cm"'') (i.e. above which most material is in protostars). The dif- 
ferent lines are as in Figure[3] Thick lines denote those calculations with ra- 
diative feedback, while thin lines are using the barotropic equation of state. 
The line types and colours denote the magnetic field strength (also ordered 
from top to bottom in each panel with progressively increasing magnetic 
field strength). Magnetic fields can be seen to affect the collapse rate at all 
density thresholds (all panels), while radiative feedback primarily prevents 
fragmentation in the highest density regions of the cloud (top panel, com- 
paring thin and thick lines). 

though the overall masses are lower) apart from some divergence at 
t > lAta in the M/$ = 3 calculation. By contrast, at higher den- 
sities (p > 10~^'^g cm~^, top panel), where the gas is optically- 
thick to radiation, there are differences of up to ~ 50% in the mass 
above this density between the barotropic and radiative transfer cal- 
culations (the latter having systematically lower mass accumulation 
rates) similar to the differences observed in Figure[3] 

4.4 Dynamics 

The effect of the reduced fragmentation on the dynamics of the 
protostars due to the radiative feedback is illustrated in Figure [5] 
showing a time sequence of the evolution in the strongest magnetic 
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Figure 5. Time sequence of fragmentation in tlie strongest magnetic field calculation (Af/<1? = 3), shown from the comparison time of i = 1.2 iff in Figure|2] 
(top row, where the panels here are shifted in position relative to Figure|2]to follow the subsequent fragmentation) up to t = 1.45fj (bottom row) at intervals 
of O.OStff . The two left-hand columns show column density for the calculation using a barotropic equation of state (first column) and with radiative transfer 
(second column). Corresponding mass-weighted temperatures for the two calculations are shown in the two right-hand columns, highlighting the heating of 
the gas due to the radiative feedback in the regions immediately surrounding the protostars. Although the radiative feedback suppresses fragmentation on the 
smallest scales, dynamical interactions nevertheless occur over larger length scales, as evident from the merger of the two star-forming cores which occurs at 
t = lAtff in each of the calculations. 
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field case (M/$ = 3) in intervals of 0.05 ts (9,500 yrs) from the 
onset of star formation. As in Figure[2l the two left-hand columns 
show column densities for the barotropic (first column) and radia- 
tive transfer (second column) simulations, whilst the corresponding 
right-hand panels show the mass-weighted temperature. The over- 
all picture is similar to that apparent from Figure |2l except that the 
time sequence shows that despite the fact that the Af/"I> — 3 sim- 
ulations each only produced a single object at f = 1.2 collapse 
and fragmentation continues to form binary and multiple systems 
at later times (i.e. t > 1.3 1//). However, whilst the larger scale dy- 
namics is similar between the barotropic and RT calculations (e.g. 
the merger of the two main collapsed regions at f ~ 1.4 tg), the 
smaller scale disc fragmentation which results in a triple ejection 
at t = l.Stff in the barotropic calculation is completely absent 
from the radiative transfer simulation, as a result of the higher tem- 
peratures surrounding the first protostars to form (fourth column of 
Figure O. 

The combined effects of stronger magnetic fields and radiative 
feedback are, therefore, to decrease the star formation rate (Figure 
O and decrease the number of occurrences of dynamical interac- 
tions and ejections between protostars. 

4.5 Protostellar masses 

The calculations presented here are of 50 M0 clouds that collapse 
to form 3-23 protostars. With such small numbers of objects, and 
the fact that the simulations have not all been followed for the same 
amount of time, it does not make sense to attempt to plot stellar 
mass functions. Rather, in Table[T] we give the amount of mass that 
has been converted into protostars (sink particles), the number of 
protostars, and the mean and median masses of the protostars. We 
give these values at f = 1.25 ta for all but one of the calculations, 
and at the end of each calculat ion. 

Generally, as found by ISatd ( l2009bh . the effect of radia- 
tive feedback is to dramatically decrease the number of proto- 
stars formed compared with the barotropic equation of state (Table 
[H columns 3 and 8). Simultaneously, the protostars are generally 
found to be more massive with radiative feedback because gas that 
would have formed other objects via the fragmentation of discs and 
nearby filaments using a barotropic equation of state is hotter and 
is able to be accreted by existing protostars instead (Table [T] col- 
umn 5 for A//$ = 5, 10 and column 8 for M/^ = 00). In the 
strongest magnetic field case, these statements are still true, but the 
trends only become apparent fairly late in the calculations because 
of the delay of the star formation caused by the strong field (Table 
□ column 8 for M/$ = 3). 

When investigating the effect of the magnetic field things be- 
come more interesting. As already discussed, the rate at which gas 
is converted into stars decreases strongly with increasing magnetic 
field strength for both the barotropic and radiative transfer calcula- 
tions (Table [U column 4). However, where this mass goes differs 
significantly between the barotropic and radiative transfer calcula- 
tions. For the barotropic calculations, the rate of protostar forma- 
tion decreases strongly with increasing magnetic field strength, but 
the typical masses of objects are independent of the magnetic field 
strength (Table[T] columns 5,6, 10, and 11). We also note that if the 
calculations are followed for a long periods of time all of the cal- 
culations eventually produce large numbers of objects regardless 
of the field strength (Table [T] column 8). However, with radiative 
feedback there is no significant dependence of the rate of protostar 
formation on the magnetic field strength (Table [T] columns 3 and 
8), and there is an indication that the mean masses of the proto- 



stars may increase with decreasing magnetic field strength (Table 
[T] columns 5 and 10). This latter effect is presumably because more 
of the gas is supported with a stronger magnetic field and not able to 
be accreted by the protostars. Although this needs to be confirmed 
with larger calculations that form larger numbers of objects, this 
implies that the characteristic stellar mass may decrease with in- 
creasing magnetic field strength, a result that is somewhat counter- 
intuitive since a naive calculation of a magnetic Jeans mass would 
lead one to conclude that the characteristic stellar mass should in- 
crease with increasing magnetic field strength. 



5 DISCUSSION 

In this paper we have studied, for the first time, the combined ef- 
fects of magnetic fields and radiative feedback on the formation of 
stellar clusters from turbulent molecular clouds. We find that the 
two effects are complementary in the sense that they affect the star 
formation process at very different scales. Magnetic fields affect 
the large-scale cloud structure (Figure [T]l, influencing all densities 
in the cloud (Figure lU. Stronger fields decrease the overall star 
formation rate (Figure [Sj. By contrast, radiative feedback affects 
small-scale fragmentation (Figure[2l( and influences only the high- 
est densities in the cloud (Figure©. It influences the star formation 
rate primarily by inhibiting small-scale fragmentation in cores once 
the first protostar has been formed (Figures (2] [3] and [5j. However, 
multiple systems are still common, formed from well-separated but 
mutually bound condensations (Figure[5j. 

The primary effect of the magnetic field is to lower the accre- 
tion rate onto the star-forming cores by providing large-scale sup- 
port to low-density regions of the cloud, thus preventing this ma- 
terial from subsequently being accreted. There is no clear shift in 
the onset of star formation with magnetic field strength (Figure[3](, 
except perhaps in the strongest magnetic calculation where star for- 
mation (ie. sink particle creation) does not initiate until t ~ 1.20tft 
(Figures [3] and [Sj compared to t ~ 1.03 — l.lltfl in the moder- 
ate/weak/zero field simulations. Since the simulations do not pro- 
duce large numbers of protostars and they are not followed very 
far any conclusions regarding the masses of the protostars must 
be treated with caution. However, we find that using a barotropic 
equation of state the typical masses of the protostars do not depend 
significantly on the magnetic field strength but the number of pro- 
tostars formed increases with weaker fields (Table [TJ. Conversely, 
with radiative feedback, the numbers of protostars formed in the 
clouds does not vary greatly with the magnetic field strength but 
the masses of the protostars tend to be lower with stronger magnetic 
fields. Generally, radiative feedback results in a larger characteristic 
protostellar mass than using the barotropic approximation. 

The general effect of radiation on the fragmentation is eas- 
ily understood in terms of the increase in the Jeans length of the 
heated gas surrounding existing protostars. An increase in temper- 
ature (e.g. from lOA" to > 30A' as in Figure|2l( leads to an increase 
in the Jeans length since Xj oc T^^'^. Because radiative feedback 
acts mainly on small scales, it takes longer for this to affect the 
overall star formation rate substantially (e.g. note the reduction in 
the figures of M{> p) propagating slowly to higher densities in 
the Af/"I> — 3 run in Figure |4] due to the progressive heating of 
wider regions of the cloud visible in Figure [3} . However, it has a 
dramatic influence on the initial mass function (IMF) by suppress- 
ing fragmentation in discs (and nearby filaments) and decreasing 
the likelihood of forming multiple systems from which low-mass 
members can be ejected (e.g. Figure [Sjl. The effect of the radia- 
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0.06 




5 
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0.11 


0.12 
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0.22 


0.22 


0.22 


1.525 


18 


1.96 


0.11 
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1.235 


10 


2.09 


0.70 


0.78 




10 


2 


0.75 


0.38 


0.38 


1.362 


5 


1.92 


0.38 


0.14 




5 


4 


0.50 


0.13 


0.13 


1.437 


10 


2.34 


0.23 


0.21 




3 


1 


0.21 


0.21 


0.21 


1.541 


7 


1.80 


0.26 


0.18 



Table 1. The statistical properties of the protostars formed in the eight calculations. For each of the four mass-to-flux ratios, barotropic and radiative transfer 
calculations were performed. Due to computational expense, the calculations were followed for a different amounts of time. All but one calculation was 
evolved until 1.25 tg, so we give the statistical properties of the simulations at this time. We also give the statistical properties at the end of each calculation. 
In each case, we give the number of protostars (sink particles) formed, the total mass in protostars, and the mean and median masses of the protostars. It is 
clear that using a barotropic equation of state produces many more objects than are obtained with radiative transfer It is also clear that the rate of protostar 
production decreases strongly with magnetic field strength in the barotropic calculations (column 3), although the typical mass of the protostars is independent 
of field strength. Conversely, with radiative transfer there is no significant dependence of the rate of protostar production with magnetic field strength, and 
there is an indication that the mean masses of the protostars may decrease with decreasing field strength (columns 5 and 10). 



tive feedback is also more pronounced when the potential well in 
which the protostars form is deeper, partially offsetting the inherent 
decrease of the Jeans mass with increasing density, which in turn 
leads to a reduced dependence of the IMF on th e initial density of 
the cloud which, as discussed in detail by iBatd ( l2009bh . may ex- 
plain why the IMF appears to be so universal across very different 
star-forming environments. 



5.1 Comparison with observations 

With a cloud of only 50 solar masses it is difficult to make a 
statistically meaningful comparison with observed star-forming 
molecular clouds as a whole since nearby clouds typi cally contain 
10^- 10'' M0 of material over areas as large as 74 pc^ tevans et aP 
12009 *). Rather our simulated clouds fall within the definition of a 
'millimeter core', ie. uyi^ >^ 2 x 10* cm^"^, and sizes similar to 
the typically measured core sizes of 1.5 x 10"^ AU in Ophiu chus 
and 3 x IC^A U in Perseus and Serpens tenoch etalJboOTh . As 
pointed out bv lBateetalH2003h . the dense cores formed in sim- 
ulations of the size presented here are similar to the Ophiuchus-F 
core which measures ~ 0.1 pc across and has a mass of 8 Mq 
jMotte, Andre & NerillT998h . In Figure |6] we show simulated ex- 
tinction maps of the four runs with radiative feedback , at the sim- 
ulation resolution (left) and at the resolution of the lEvans et"aD 
( l2009l) extinction maps for Ophiuchus (right), with extinction on 
a linear greyscale map fro m Ay = 1 t o Ay — 25 which may 
be directly compared to the lEvans et al ] l l2009h maps. To produce 
the extinction maps we have simply used the inverse of the con- 
version from extinction to hy drogen column de nsity of 1.37 x 
10^^cm~^mag^^ adopted bv lEvans et alj l l2009l) . The resolution 
of the right hand panels is calculated from the c2d extinction map 
resolution of 270", which at the assumed distance of 125pc for 
Ophiuchus gives a resolution of 0.16 pc. We have simulated this 
resolution in our maps by enforcing a minimum smoothing length 
of O.OSpc on the SPH particles when calculating the column den- 
sity (i.e., approximating the point spread function for the extinction 

maps by the SPH kernel smoothing function). 

Star formation efficiencies are calculated by lEvans et"aD 

( I2OO9.) by dividing the mass in Young Stellar Objects (YSOs, de- 
fined as objects with infrared excesses assumed to correspond to 



the presence of a disc) by the total mass of the cloud plus YSOs. 
That is, 



SFE : 



M, + M{doud) ' 



(21) 



where M{cloud) is derived by integrating the extinction maps, 
converted to column density, over area. Despite the low resolu- 
tion of the observations compared to our simulated cloud, the cloud 
masses measured from the clouds on the right hand side of Figure[6] 
by integrating column density over the area within the Ay = 2 
contour, are remarkably accurate. For example, the measured mass 
for the zero magnetic fields case (Af/"I> = 00) at the observa- 
tional resolution is 44. 1 Mq , which may be compared with the total 
cloud mass in our simulations of 50 M©, of which ^ 44 Mq lies 
within the sphere with the approximate radius of the Ay = 2 con- 
tour. The caveat to this for the observations is that the conversion 
from extinction to column density relies on a model for the dust, 
chan ges to which can have a significant impact on cloud masses 
(e.g. lEvans etal]|2009l discuss the fact that their cloud masses are 
revised down by a factor of 1.4 from previous estimates due to re- 
vision of the dust model). 

Efficiencies thus derived by lEvans et"ah I 12009) range from 
3 — 6%, which is assumed to represent an average over the last 
2 Myr given that this is the estimated lifetime of YSOs with infrared 
excesses. The comparison between observations and our simula- 
tions is made more difficult by the fact that we are not able to 
follow the cloud collapse for longer than around 1.5tff with cur- 
rent computational resource^ Nevertheless, one can make tenta- 
tive estimates based on the star formation rates we find in Figure [3] 
and the masses in Table [T] For the four clouds shown on the right 
hand side of Figure |6| (the four runs with radiative feedback in- 
cluded), cloud masses measured by integrating the column density 
within the Ay = 2 contour are 44.1, 44.0, 43.6 and 43.9 for the 
M/^ — 00, 10, 5 and 3 clouds respectively. A straightforward 
application of i21\ at the end point of each of our calculations indi- 
cates that of order 3 — 5% of the gas has been converted into stars 

^ The key limitation being that for only a few collapsed objects, good load 
balancing is very difficult to achieve, limiting the usefulness of simply run- 
ning on a higher number of processors. 
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Figure 6. Extinction maps of the four cal culations wi th radiative transfer and varying initial magnetic field strength (A//<i> as indicated) at the simulation 
resolution (left) and at the resolution of the lEvans et alj J200 9) extinction maps for Ophiuchus (right). The Ay = 2 contour is shown, the mass inside of which 
was used to calculat e the mass in the c2 d survey. We have converted hydrogen column density to Ay using the inverse of the conversion factor adopted for the 
c2d observations bv lEvans et J2009l) . i.e. a conversion from extinction to hydrogen column density of 1.37 X lO^^cm ^mag ^. The scale from Ay = 1 
to Ay = 25 is the same as that used for the c2d maps. 



over the time for which the simulations have been run. Whilst these 
values are in agreement with the observational results (though not 
for dense gas), they are not very meaningful given that they repre- 
sent evolution over fractions of a freefall time beyond initial star 
formation (the end time for each of the calculations is given in Ta- 
ble[T]and can be inferred from FigureO and will increase with time 
as more mass is converted into stars. 

More useful are estimates which take into account the 
timescale over which star formation has proceeded. The depletion 
time for the cloud is given by 

tdep = M (cloud) /M, . (22) 

Calculating the average star formation rate from the onset of star 
formation using Figure [3] and assuming that these rates will con- 
tinue indefinitely (a dubious assumption), we obtain depletion 
times of 0.8, 1.6, 1.3 and 1.8 Myr for the clouds in the above 
four calculations respectively. Whilst these are very short com- 
pared to the global depletion times for the clouds in 'Evans et al.' 
(2009) of 30 to 66 Myr, they are in agreement with the depletion 
timescale derived for dense cores within such clouds (i.e., gas with 
n > 2 X lO^cm^'', which our initial cloud density lies above) 
which are in the range of 0.6-2. 9 Myr with an average of l.SMyr. 

Finally, lEvans et al ]( l2009l) quantify the observed inefficiency 
in terms of the star for mation rate per free-fall time, defined as 
jKrumholz & Tanll2007h 

SFRs = M,ts /Mcioud, (23) 

where tg is defined as the free-fall time for the mean density of the 
cloud and which here we take as the initial free-fall time for our 
initially uniform density clouds. 

Using this measure we find star formation rates of SFRs = 
0.23, 0.12, 0.15 and 0.10 for the four runs (M/<E> = oo, 10, 5 and 
3 respectively) that include radiative feedback and SFRa — 0.32, 
0.18, 0.17 and 0.12 for the four runs using a barotropic equation of 
state, that is, neglecting radiative feedback. Thus, only the strong 



magnetic field calculations (/? < 1, corresponding to Af/<I> < 7) 
that include radiative feedback ap proach the observe d range of 
SFRff = 0.03 - 0.06 measured by fevans et all ( l2009h . All of the 
calculations with weaker field strengths and/or neglecting radiative 
feedback have star formation rates that are much higher than ob- 
servations suggest. From the point of view of matc hing theory to 
observ ation, this is reassuring, since, as discussed in lPrice &Batel 
the most realistic of our calculations in terms of magnetic 
field strength is the strongest field case, — 3, since molecu- 

lar cloud cores are typically observed with mass-to-flux ratios that 
are marginally supercritical (i.e. M/$ ~ 2 — 3) and with magnetic 
pressure smaller than gas pressure by a f actor of ~ 3 (i.e., /3 ~ 0.3) 
(ICrutcheill 19991 : iHeiles & TrolMi3l2004) . However we caution that 
any conclusions regarding the star formation efficiency from these 
calculations are necessarily limited by the relatively short period 
over which we have been able to follow the calculations beyond 
the free-fall time. 

Furthermore our results present only a lower limit on the effect 
of feedback since we have neglected feedback from within 0.5 AU 
of a star including the driving of stellar winds and collimated out- 
flows which may act to fii rther reduce the star formation efficiency 
jMatzner & McKedl200ol) . perhaps explaining the remaining dis- 
crepancy between the efficiencies we find and the observed range 
of 3 -6%. 



5.2 Implications for theory 

The reduction in star formation rate is primarily a result of the sup- 
port provided to the cloud by the magnetic field. The global mag- 
netic field, whilst not sufficient to prevent collapse altogether, is 
neve rtheless able to aff ect the binding energy. 

ICIark et al.l ( I2OO8 ) point out that the star formation rate can 
be made arbitrarily low in globally unbound clouds by increasing 
the ratio of kinetic to gravitational potential energy Ei:in/\Egrav\ 
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(set to unity in the initial conditions for the calculations we present 
here). The fact that increasing the turbulent velocity dispersion 
can decrease the efficiency of star formation in the sense of low- 
erin g the star formation rate has also been discussed previously 
(e.g JPadoadl 19951 : llflessen et aLlllOOOh . However, the kinetic en- 
ergy cannot be increased indefinitely for a cloud of this size with- 
out violating the observational constraints on the turbulent veloc- 
ity field. Observationally, the velocit y line w idth scales with cloud 
size a pproximately as v o c dLarson 
1 19871 : lBrunt&Heve3l2002l : iHever & Brunt 



Isolomon et aP 
2004) with a magni- 



1981 



tude of V ~ 1 km s on 1 pc scales and a scatter of a factor 
of two ( iHever & Brun3l2004h . For a cloud the size of those mod- 
elled here (0.375 pc), this gives a typical velocity dispersion of 
V ~ 0.6 km s~^ (Mach 3.3) which is almost a factor of two less 
than the velocity dispersion of our initial conditions. Thus, our ini- 
tial conditions are already at the upper end of the observed velocity 
dispersion in molecular clouds so there would appear to be little 
scope for achieving a lower star formation rate by boosting the level 
of turbulence. 

By contrast, as we have shown through the simulations pre- 
sented here, a low star formation rate requires only a magnetic field 
of similar st rength to observ ational estimates (i.e. a mass-to-flux 
ratio of > 3, ICrutcher|[T999h and the effects of radiative feedback 
which has no large free parameters (once the metallicity is set). 
Similar results with regards to the reduct ion in star formation rate 
with m agnetic field strength are found bv lVazquez-Semadeni et all 
( I2OO5 I) in the context of (scale-free) driven turbulence simulations. 

It therefore appears that both strong magnetic fields and radia- 
tive feedback from protostars are crucial ingredients in regulating 
star formation to a slow and inefficient level, which cannot be ne- 
glected from numerical simulations of the star formation process. 
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